% 火箭公式计算
% 参数是随便写的，不具有实际意义
% Gitee Repo

clc
clear

m0 = 10; %初质量
mfinal = 2; %末质量
v0 = 0; %初始速度
u =  20; %燃料速度

dt = 0.001; % 时间步长

dmdt = 10; %燃料喷射速度
dm = dmdt*dt; %每次喷出的燃料质量
g = 10; %重力加速度
f = 1; %摩擦力系数

N = ceil((m0-mfinal)/dm);
t = zeros(N,1);
v = zeros(N,1); %各时刻速度，考虑重力和摩擦力
videal = zeros(N,1); %各时刻速度
h = zeros(N,1); %各时刻高度，考虑重力和摩擦力
hideal = zeros(N,1); %各时刻高度
m = zeros(N,1); %各时刻质量

v(1)=v0;
videal(1)=v0;
h(1)=0;
hideal(1)=0;
m(1)=m0;
t(1)=0;


for tick = 2:N
    t(tick) = t(tick-1)+dt;
    m(tick) = m(tick-1)-dm;
    dv = u*dm/(m(tick)); %速度增量

    videal(tick) = videal(tick-1) + dv;
    v(tick) = v(tick-1) + dv - g*dt - f*v(tick-1)*dt;

    h(tick) = h(tick-1) + v(tick)*dt;
    hideal(tick) = hideal(tick-1) + videal(tick)*dt;
end
disp(videal(N))
disp(u*log(m0/mfinal) )
disp(v(N))

hformula = u*log(m0./(m0-dmdt*t));

figure
hold on
##axis equal
plot(t,videal,t,v);
plot(t,hformula,'-');
legend('v ideal','v with gravity and friction','v formula')

